8 Zooarchaeology
8.1 Case studies
The following map shows the sites under investigation, divided by chronology. Please select the desired chronology (or chronologies) from the legend on the right. Legend: R = Roman, LR = Late Roman, EMA = Early Middle Ages, Ma = 11th c. onwards
The faunal dataset used in this study is both extensive and diverse, containing over 466 records. While NISP is a useful proxy for historical animal farming, it is not without its limitations, as previously discussed in the methods section. What is crucial to emphasize here, however, is the presence of overdispersion within the data, which necessitates more sophisticated and nuanced approaches than simply calculating overall means for each animal species. Overdispersion is a common occurrence when analyzing datasets of this nature, given the unique factors at play in each context, including historical and depositional influences. Nevertheless, data modeling requires simplification and causal reasoning, and the best approach is to start with straightforward models that can account for overdispersion and generate credible distributions. To this end, Bayesian hierarchical models were developed for each chronology, context type, macroregion, and geography. As a further step, an additional analysis was conducted that solely focused on altitude and chronology as predictors. By examining the specific contributions of altitude and chronology in predicting animal farming/consumption patterns, this analysis provides valuable insights that complement the earlier models. These findings underscore the importance of considering multiple factors when studying the probability of occurrence of farmed and wild animals in historical contexts. They also demonstrate the potential benefits of simplified models to focus on key predictors. Moreover, several attempts were made to create a coherent understanding of the likelihood of economically valuable animals occurring during the first millennium, in order to provide a more comprehensive perspective on the role of animal farming in shaping historical societies.
8.2 Data exploration
As previously mentioned, the faunal dataset used in this study exhibits overdispersion, which is a common issue when analysing datasets of this type. The term ‘overdispersion’ refers to the presence of greater variability in the data than expected based on a normal curve (Figure 8.1). In this section, I will present the distribution of the animals of interest to provide a visual representation of the dispersion within the dataset. By examining the distribution curves, we can better understand the variability that exist within the dataset, and can use this information to develop more accurate models. In Figure 8.2, it is evident that the distribution of animal remains in the faunal dataset is not symmetrical and does not conform to a normal curve. This non-normal distribution indicates that the standard measures of central tendency, such as mean and median, may not accurately capture the overall distribution of the data. Traditional statistical measures, such as measures of central tendency and dispersion, are commonly used in frequentist approaches to analyse data. However, these measures may not be appropriate for analysing the complex patterns of animal farming and consumption in the faunal dataset due to the presence of overdispersion. As an alternative, Bayesian multilevel models can account for overdispersion by incorporating appropriate probability distributions, such as the betabinomial distribution.


8.3 Chronology
8.3.1 Trends by century
The proposed model uses a betabinomial distribution to model overdispersion in the data and estimates the precision (or shape) parameter \(\phi\) in the Beta distribution. The \(A\) on the left side of the formula is the outcome variable, which represents the animal NISP count for each observed sample \(i\). The model is intercept-only, meaning that there are no further predictors included in the model. The average probability of success \(\bar{p}_{i}\) is modelled through a logit link function, and is assumed to equal the intercept \(\alpha\). The model is also indexed with \({[Century]}\) so that separate estimates can be obtained for each century (1st BCE to 11th CE).
\[ A_{i} \sim BetaBinomial(NISP_{i}, \bar{p}_{i} , \phi_{i}) \]
\[ logit(\bar{p}_{i}) = \alpha_{[Century]} \]
\[ \alpha_{[Century]} \sim Normal(0,1.5) \]
\[ \phi_{[Century]} \sim Exponential(1) + 2 \] The prior distributions selected for the model are weakly informative. Specifically, the prior chosen for the intercept \(\alpha\) is a normal distribution with a mean of 0 and a standard deviation of 1.5. Given that \(logit(\bar{p}) = \alpha\) is modelled using the logit link, a mean of 0 corresponds to a probability of 0.50 on the logit scale. With a standard deviation of 1.5, the prior is flat and only slightly informative.
As for the shape parameter \(\phi\), the prior was selected to ensure that the parameter was at least 2, which corresponds to a flat distribution. This value (2) was then added to an exponential distribution with rate 1 to produce a weakly informative prior.
A prior predictive simulation (Figure 8.3) shows how when no information is provided to the model, counts on the extremes of the binomial histogram are more likely, while all other counts have similar probabilities. Bimodal priors are common practice with Beta distributions as they are non-informative enough. However, updating a bimodal prior even with one observation rapidly changes the posterior.


The model for pigs’ NISP counts shows a positive trend above the millennium average (0.36) for predicted high density intervals (HDIs) from the 1st century BCE to the 3rd century CE. The probability intervals exhibit a steady decrease that begins in the 3rd century and continues until the 8th century, when the mean of the probability distribution coincides with the first millennium average. Despite positive trends reported in the 8th and 9th centuries, interpretation of these values requires caution since the credible intervals are wider due to fewer observed samples dated to these centuries. The HDIs mean for the 10th and 11th centuries is again around the millennium average, indicating a decrease in the probability of pig occurrence.
Cattle, on the other hand, exhibit different trends, with a HDI interval in the 1st century BCE above the millennium average (0.20), although with a large credible interval, and a negative trend from the 1st to the 3rd century CE. HDIs are again increasingly positive between the 4th and 7th centuries, but decrease in the 8th century. As for pigs, HDIs means are close to the millennium average in the 10th and 11th centuries.
The modelled probability of occurrence of caprine reveals a positive trend in the 1st century BCE, followed by a negative trend (below the average of 0.28) from the 1st to the 3rd century CE. The probability of occurrence then increases around the 4th/5th centuries up to the 6th/7th centuries, with HDIs stable around the mean from the 8th to the 11th century.
Edible wild mammals and domestic fowl predicted probabilities are relatively stable trends throughout the analysed period, with HDI means consistently below 0.10. Domestic fowl HDIs remain below the millennium mean of 0.06 from the 1st century BCE to the 2nd century CE, and then increase steadily from the 3rd century to remain around the mean until the 7th century. From the 7th to the 11th century, the HDIs remain positive, but with a much larger credible interval. Similarly, wild mammals’ HDIs are positive in the 1st century BCE and remain stable around the millennium average (0.05) until the 6th century CE, after which they also show a slight increase until the 11th century.
In addition to plotting the probabilities of occurrence extracted by the models, the distributions of the precision parameter \(\phi\) are plotted for samples in each century. Each line is plotted with two probability intervals 0.65 (shorter thicker line) and 0.99 (longer thinner line). The probability density functions are also added to the top of the bar.


8.3.2 Trends by phase
The relationship between animal NISP counts and chronological phases is similar to the previous one, but with a different indexing variable for the intercept. Instead of being indexed by century, the intercept \(\alpha\) is now indexed by phase (\({[ChrID]}\)). While a phase can span several centuries, this modification allows the model to estimate separate intercepts for each phase, providing more specific insight into potential differences in the outcome variable and the average probability of animals occurrences over time. The rest of the model structure and distributional assumptions remain the same as the previous model, which uses a betabinomial distribution to model overdispersion in the data and an intercept-only structure.
\[ A_{i} \sim BetaBinomial(NISP_{i}, \bar{p}_{i} , \phi_{i}) \]
\[ logit(\bar{p}_{i}) = \alpha_{[ChrID]} \]
\[ \alpha_{[ChrID]} \sim Normal(0,1.5) \]
\[ \phi_{[ChrID]} \sim Exponential(1)+2 \]
As anticipated, the phase-level trends exhibit posterior distributions similar to the more detailed century-level results. However, the phase-level model offers an advantage - the credible intervals are considerably narrower than those of the previous model. Since each phase spans several centuries, except for the medieval phase that only includes the 11th century, there are more data points available to provide more reliable probabilities. Pigs’ NISP are at their highest point during the Roman phase, whereas they decrease in the Late Roman and Early Medieval phase, only to be slightly below the millennium average (of 0.37) in the 11th century. The predicted millennium average is comparable to the previous century-level model (0.36), with a small difference which can be due to the grouping of the observations. The diachronical means for the other animals are also very comparable to the other with a maximum difference of ±1 for each animal. Cattle’s NISP probabilities of occurrence increase throughout the millennium, although in the 11th century (marked as Ma on the plot) the 99% credible interval is larger, spanning from 0.15 to 0.26. Conversely, the century-level model for this period indicated a small increase towards the millennium average. As for the century-level model, the caprine are the second most represented animal in the millennium. Their probability of occurrence shows HDIs just below the millennium mean (of 0.27) in the Roman and Late Roman age, and positive HDIs values for the Early Medieval and Medieval phase, a situation not so different than the century-level model. The phase-level trends for domestic fowl and wild mammals show perhaps clearer trends than the century-level models. The HDI for domestic fowl, which is below the millennium mean of 0.05 in the Roman age, slightly increases in the later phases, with a peak in the 11th century. This peak is however less reliable than other phases, as the credible interval range is much larger. Wild mammals posterior probabilities on the other hand are stable around the mean (0.04) with a small increase in probability during the Early Medieval and Medieval phases. In addition to improved credible intervals as opposed to the century-level model, the probability intervals for the \(\phi\) shape parameter are also more reliable, as the PDFs are more narrow and most of the values are around the mean. Unfortunately, a smaller amount of samples for the 11th century produced more dispersed curves.


8.3.2.1 Community plot
After estimating the intercepts for each animal separately, they were compiled together to get a cohesive view of the probability of occurrence of the four animals. This approach is not ideal, but it can still be informative with regards to the patterns of animal exploitation in the archaeological record. The joint community plot allows to visually compare the probability of occurrence of each animal species over time and identify trends in their occurrence. For instance, the plot below shows how the predicted probability values for the three main domesticates—pigs, cattle and caprine vary across the four phases. Pigs have higher probabilities during the Roman phase, while already in the Late Roman period they get closer to other domesticates. It is in this period that the inter-sample variability of \(\phi\) seem to decrease, although the values are not very high. The precision \(\phi\) is instead very variable for edible wild mammals and domestic fowl, as it can be expected given that it is much rarer to find these animals in a sample for various reasons already discussed in the methods. The high variability in the precision of the medieval (11th century) intercept and \(\phi\) estimates can be explained by the lower number of samples. Overall, the most noticeable trends in the plot are the increased number of cattle in the early medieval period and the decreased number of pigs from the Late Roman phase. Modelling each animal species separately does not fully capture the complex interactions that exist among them, and future work should explore more advanced modelling approaches that account for the co-occurrence of multiple animal species in the same samples.

8.4 Context type
Chronological models provide a broad overview of the trends in domestic and wild animal husbandry and consumption. However, it is important to consider the context where the animal remains are found because different site types can influence the NISP of a particular animal. We can stratify the dataset by site type, but there is a challenge: the four chronologies we are studying also affect the availability, frequency, and economic dynamics of the site types. For example, the agricultural strategies employed by rural villas during the Roman and Late Roman periods could differ based on political and economic events, which could affect the relative amount of animal NISP. To better understand how the context where animal remains were found affects our outcome variable, the NISP, we can use a Directed Acyclic Graph (DAG). In this DAG, the Chronology node is a parent of the Site_Type node, and it affects the NISP count directly and indirectly through Site_Type. To account for this backdoor path, we block by Chronology and stratify the dataset by Site_Type. Since we have two categorical predictors of interest, we can use an interaction index variable to generate distinct intercepts for each chronology and context type. For instance, the Roman:Rural context type could be assigned an index of 1, while Roman:Urban could be assigned an index of 2, and so on. This interaction dummy index (\({[TCid]}\)) will determine the variation in intercepts (\(\alpha\)) across different context types.

The approach for the proposed model for estimating the posterior likelihood of animal occurrence in various chronological phases and context types is very to the models discussed earlier. The only modification is the interaction index. The model structure and priors remain the same as in the previous models. There are several contexts that have been recorded in the database (see Table 5.1 for a full list), but the model required some simplifications in order to avoid extremely small categories. For instance, religious monasteries have been grouped with other Religious sites, or Roman mansiones have been classified as Rural for these models. Categorical simplifications can obscure the nuances of data, but tiny sample sizes produce very uninformative posterior distributions.
\[ A_{i} \sim BetaBinomial(NISP_{i}, \bar{p}_{i} , \phi_{i}) \]
\[ logit(\bar{p}_{i}) = \alpha_{[TCid]} \]
\[ \alpha_{[TCid]} \sim Normal(0,1.5) \]
\[ \phi_{[TCid]} \sim Exponential(1) + 2 \]
8.4.1 Pigs
As previously noted, pig remains are consistently the most frequently found type of faunal remains in first millennium Italian excavations. However, this study reveals divergent patterns in the presence of pigs depending on the site type and function. We will first discuss categories with narrower credible intervals before moving on to those with smaller sample sizes and wider credible intervals. Urban sites are the most common category, where the probability of finding pig remains follows patterns similar to the chronological trends presented earlier. During the Roman phase, the probability of finding pig remains in urban contexts is very high compared to other animals or categories, with a mean probability of almost 0.5. The consumption of pork in cities appears to have decreased during the Late Roman and Early Medieval phase, with similar credible intervals, only to increase again in the 11th century. Possible hypotheses will be discussed later on. Although the precision of the pig beta posterior in urban contexts is slightly below the across-contexts mean of 5, the curves are not too dispersed, at least in the Roman and Late Roman periods, and we can be confident in these results. On the other hand, the 95% HDIs for pigs in rural contexts are lower than those in urban contexts during the Roman and Late Roman periods. In the Early Medieval phase, the mean probabilities are comparable to urban contexts, with both categories having a mean of around 0.31, which then increases to 0.37 in the 11th century (0.39 in urban contexts). In general, while the decline in pig presence appears to be significant in urban areas, rural contexts exhibit only a slight increase. Rural villas have been analyzed separately from other rural sites due to their potential implications for the production and consumption of meat, given their status as elite sites. Although the sample size is not large enough to provide highly confident posterior predictions, the results suggest that pork was likely an important component of the diet in these sites, with a mean of approximately 0.4 in the Roman period and an increase in the Late Roman phase. However, the Early Medieval phase is characterized by greater uncertainty, as the sample size is reduced and rural estates often changed function or were abandoned. The credible interval for this period ranges from 0.25 to 0.62. The suggested increase in pigs NISP in the Late Roman period will be discussed further later, as the historical debate might provide helpful hypotheses. It is worth noting that there is no sample from the 11th century for rural villas, which is why that line is not included in the graph. Fortified sites were also large consumers of pigs, with the probabilities of occurrence sticking around the mean value and a minor increase in the Early Medieval-Medieval periods. However, there are no significant changes in trends to observe. The probabilities of occurrence in the Late Roman phase are not so trustworthy as the credible interval is large. It is worth noting that there is no Roman posterior prediction on the graph because there were not fortified Roman sites in the sample except for one, the castrum of Ostia (MacKinnon, 2014), which has not been included in the observed data for this model as it was one single sample. The castrum had a 71.9% NISP proportion of pigs, out of a sample of 121 total NISP. The religious category has fewer observations than the site types presented before, and the credible intervals are wider. The Roman and Late Roman religious sites mostly include temples such as the temples C and D of Grumentum, the Demetra temple in Macchia delle Valli, the Mithraeum of Crypta Balbi, and more. In the Medieval phases, religious sites are mostly monasteries including the well-studied ones of Monte Gelato, Farfa, San Salvo, S. Giulia of Brescia, and San Vincenzo al Volturno. The range of credible intervals can span probabilities of 0.40, so we must be cautious in our considerations. If we look at the mean of the HDIs, a precise pattern does not emerge. However, in every chronology, the probability of occurrence of pig’s NISP is over the across-contexts average of 0.41. The credible intervals for necropolis sites are too large - spanning almost the entire probability range - to trace any clear patterns. With only nine samples from seven sites, including the Roman and Late Roman necropoleis of Cantone, Otranto, San Cassiano (Riva del Garda), Trieste (loc. Crosada), Poggio Gramignano, San Lorenzo di Sebato, and the Early Medieval necropolis of Baggiovara, the sample size is small. As a result, it is best to draw only qualitative conclusions in the later discussion. Finally, the last category includes only one well-known site, the Flavium amphitheater in Rome, which provided 29 zooarchaeological samples. The 95% HDIs for this site show very high probabilities, with mean values of 0.65 for the Roman and Late Roman periods and 0.49 for the Early Medieval period. However, the credible interval is quite wide for the Early Medieval phase, making it less certain. This type of site was not grouped with the other urban observations as it was too unique.


8.4.2 Cattle


8.4.3 Caprine


8.4.4 Edible W. Mammals


8.4.5 Community plot

8.5 Macroregion
To estimate the animals’ occurrences probability in each chronology and macroregion, I used a betabinomial distribution to model overdispersion in the data. The \(A\) on the left side of the formula is the outcome variable—the animal NISP counts for each observation \(i\). This is an intercept-only model, where the intercept \(\alpha\) carries an interaction index \({[REGid]}\) as the model will provide estimates for each macroregion and chronology under examination. The \(\phi\) parameter indicates the precision in the Beta distribution, modelled by chronology and macroregion.
\[ A_{i} \sim BetaBinomial(NISP_{i}, \bar{p}_{i} , \phi_{i}) \]
\[ logit(\bar{p}_{i}) = \alpha_{[REGid]} \]
\[ \alpha_{[REGid]} \sim Normal(0,1.5) \]
\[ \phi_{[REGid]} \sim Exponential(1)+2 \]
8.5.1 Pigs


8.5.2 Cattle


8.5.3 Caprine


8.5.4 Edible W. Mammals


8.6 Geography
Animals distributions can vary across different geographical features. This research has considered plain, coast, hill and mountains as the most common geographical features in the Italian peninsula. Archaeological excavations where zooarchaeological remains have been analysed are located at low altitudes. Although as expected there are more mountain sites in Northern Italy, sampled sites are evenly placed on plains, coastlands and hills across the three Italian macroregions.

To estimate the animals’ occurrences probability in each chronology and geography, I used a betabinomial distribution to model overdispersion in the data. The \(A\) on the left side of the formula is the outcome variable—the animal NISP counts for each observation \(i\). This is an intercept-only model, where the intercept \(\alpha\) carries an interaction index \({[GEOid]}\) as the model will provide estimates for each geography and chronology under examination. The \(\phi\) parameter indicates the precision in the Beta distribution, modelled by chronology and geography.
\[ A_{i} \sim BetaBinomial(NISP_{i}, \bar{p}_{i} , \phi_{i}) \]
\[ logit(\bar{p}_{i}) = \alpha_{[GEOid]} \]
\[ \alpha_{[GEOid]} \sim Normal(0,1.5) \]
\[ \phi_{[GEOid]} \sim Exponential(1)+2 \]
8.6.1 Pigs


8.6.2 Cattle


8.6.3 Caprine


8.6.4 Edible W. Mammals


8.7 Altitude
The probability of occurrence of the most common faunal remains can be modelled against the elevation of sites in the four time periods under consideration. It is worth noting that the sites where the zooarchaeological remains have been found are not evenly distributed. In the Roman age, most sites investigated are located between 0 and 100 MSL, whereas after there is an increasing number of remains from sites between 100 and 400 MSL. Whether this reflects a real shift in settlement patterns is outside the aims of this study, but it might still be informative to visualise the different distribution of sites across elevations.

The proposed model to estimate the probability of occurrence as related to the altitude (the slope \(\beta\)) and chronology (\({[ChrID]}\)) uses a betabinomial distribution to model overdispersion in the data. The \(A\) on the left side of the formula is the outcome variable—the animal NISP counts for each observation \(i\). This is a simple intercept with slope model, where the intercept \(\alpha\) carries an index \({[ChrID]}\) as the model provides estimates for each chronology under examination. A single \(\phi\) parameter indicates the precision in the Beta distribution.
\[ A_{i} \sim BetaBinomial(NISP_{i}, \bar{p}_{i} , \phi_{i}) \]
\[ logit(\bar{p}_{i}) = \alpha_{[ChrID]} + \beta_{[ChrID]}\cdot Alt_{i} \]
\[ \alpha_{ChrID} \sim Normal(0,1.5) \]
\[ \beta_{ChrID} \sim Normal(0,1.5) \]
\[ \phi \sim Exponential(1)+2 \]

8.7.1 Pigs

8.7.2 Cattle

8.7.3 Caprine

8.7.4 Edible W. Mammals

8.7.5 Community plot
